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Abstract 

Two different Reverse Monte Garlo strategies, ’RMC++’ and ’RMGPOW’, have 
been compared for determining the microscopic structure of some liquid and 
amorphous solid systems on the basis of neutron diffraction measurements. The 
hrst, 'g{r) route’, exploits the isotropic nature of liquids and calculates the 
total scattering structure factor, S(Q), via a one-dimensional Fourier transform 
of the radial distribution function. The second, called ’crystallography’ route, is 
based on the direct calculation of S{Q) in the reciprocal space from the atomic 
positions in the simulation box. We describe these two methods and apply them 
to four disordered systems of increasing complexity. The two approaches yield 
structures in good agreement to the level of two- and three body correlations; 
consequently, it has been proven that the ’crystallography route’ can also deal 
perfectly with disordered materials. This finding is important for future studies 
of liquids confined in porous media, where handling Bragg and diffuse scattering 
simultaneously is unavoidable. 

Keywords: Reverse Monte Garlo; neutron scattering; structure factor; liquids; 
modelling 


1. Introduction 

The Reverse Monte Garlo (RMC) method[i| is a simple tool used for decades 
for elucidating the detailed atomic level structure of liquids and solids from 
scattering measurements. Over the past 25 years, RMG has been successfully 
applied to a wide variety of disordered materials that display structural disorder 
of varying extent: simple liquidsH, molten salts Q, molecular Wdsd, i,0, 
waterQ and aqueous solutionsU^ metallic^ and covalent[l3, [HI glasses. A 
separate class of applications has targeted ’disordered crystals’ in which long 
range (crystalline) order and local (i.e., within the first coordination sphere) 
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disorder are present simultaneously: examples may be crystals of silver and 
copper halides [ll. [l^ and of tetrahedral molecules 141. 

It was clear early on 12| that dealing with genuine crystalline materials re¬ 
quires strategies different from those applicable for isotropic liquids/amorphous 
materials, due to the presence of long range periodic symmetries and the lo¬ 
cally anisotropic nature of crystals. Just before the turn of the millennium, the 
(so far) ultimate solution was created: the RMCPOW software 1^ is able to 
calculate Bragg- and diffuse scattering intensities directly from the particle coor¬ 
dinates, even for powder diffraction data obtained from laboratory X-ray sources 
and thermal neutron diffraction. For experimental data measured over very wide 
momentum transfer ranges, the RMCProfile strategy 16|, that involves the sepa¬ 
ration of the Bragg profile and Fourier-transform to real space, and a subsequent 
modelling of the total radial distribution function and the Bragg-profile, is also 
frequently used. The PDFGui software 17|, performing PDF-based analysis of 
powder diffraction data, is a powerful tool for providing structural models based 
on the radial distribution function of crystalline materials. This is an alterna¬ 
tive to the strictly unit-cell based investigation of crystalline structures; on the 
other hand, it is not capable of dealing with genuinely disordered structures. 
For isotropic disordered materials the original strategy of RMC[Ij may be used, 
i.e., from the atomic positions, first the radial distribution functions (RDF) are 
calculated, which later are Fourier transformed to the reciprocal space, so that 
primary experimental information, the total scattering structure factor (TSSF) 
may serve as ’target function’ of RMC. Software that can realize this strategy 
may be RMC-f-hUi, RMC_POT0 or RMC Profile [IS [ 2 ^. Details of the two 
strategies will be provided below; for now, it is important to state that a proper 
comparison between the two strategies is still missing. 

The primary aim of this work is to test these strategies for several model 
systems. Since it is obvious that the simple route, via the calculation of the 
RDF, cannot be applicable for crystals, what needs to be tested is whether the 
more time consuming ’crystallographic’ approach[l^ can be used for isotropic 
disordered systems, such as liquids. Beyond the ’per se’ interest, the timeliness 
of such a study lies in that a very important class of ’mixed’ systems, ’fluids 
in pores’ would require a method that can handle both perfect crystals (like 
zeolites) and liquids (like waterl [2li . Note that the ’crystallographic’ 

approach has already been proven to reproduce the atomic structure of sim¬ 
ple adsorbed fluids (up to the level of three body correlations) in zeolites of 
varying pore sizes using the N-RMC method in which the number of particles 
is an additional adjustable parameter!^. In that work the target structure 
factor was obtained by simulation rather than from experiments and the study 
was restricted to simple fluids. Structural investigations of such complicated 
materials, that are of utmost significance in catalysis, oil industry, soil chem¬ 
istry..., will not be possible until an established method of structural modelling 
can be proven to be applicable. Our aim now is to see whether the ’crystallo¬ 
graphic’ approach is also suitable for fitting experimental structure factors for 
more complex fluids. 

Bearing in mind the above, the two approaches are tested on disordered one 
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component systems of increasing complexity, from liquid argon to amorphous 
silicon. Liquid argon (1-Ar) is one of the simplest fluids in all respects: it can be 
easily described using radially symmetric pair potentials [ 2 ^. Liquid gallium (1- 
Ga), is a unique metallic element with possible short-lived covalent bonds that 
manifest in the slig htly unusual shape of the main peak of the total scattering 
structure factor 2^ . Liquid selenium (1-Se) is one of the most unusual elemental 
liquids, because of the twofold coordination of the atoms and the resulting 
chain-like structure[2^. Finally, amorphous silicon (a-Si) can be regarded as a 
classic example of a disordered fourfold-coordinated covalent material that, in 


contrast to its well-known crystalline form, lacks the Jon] 
the cases of 1-Ga, 1-Se and a-Si, experimental data|25 


range order [27j. In 
30l ] are from neutron 
diffraction measurements. In the case of argon, a computer-generated model of 
the liquid ( 3 ^ has been employed, for two reasons: (I) this way, no systematic 
experimental errors had to be cared for, and (2) the early experimental dataf^ 
exhibited some residual systematic errors that made a thorough comparison of 
the methods somewhat cumbersome. 


2. The two approaches for calculating the measurable total scattering 
structure factors within RMC 

Details of the RMG method can be found in various publications [ll. [Til [ 2 ^ 
[ 31 I [isl [ 1 ^ and therefore, here we will concentrate only on the parts relevant 
for calculating the structure factor from particle coordinates. 

In short, the RMC algorithm produces sets of three-dimensional particle 
coordinates for which the calculated structure factor fits the input diffraction 
data within the estimated experimental errors. The goodness-of-fit is quantified 
using a x^-value: 


Nq 




{ScalciQi') SexpiQi^^ 


( 1 ) 


where Q is the modulus of the scattering variable, the sum runs over all experi¬ 
mental points, Nq-, Sexp and Scale are the experimental and simulated structure 
factors, respectively, and cr is the ’estimated’ standard deviation for the exper¬ 
imental point i. 

To minimize j random movements are attempted for all atoms in the sim¬ 
ulation box. If the new non-overlapping position reduces differences between 
experimental and calculated structure factors, the move will be accepted. Oth¬ 
erwise, the move is accepted according to an acceptance probability, given 

by 


= min 1, exp - 


\new 


^old 


( 2 ) 
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where Xoid Xnew correspond to the original and proposed atomic coor¬ 
dinates, respectively. 

Finally, an exclusion core around each particle is defined, rcutoff, to reflect 
its effective size. If the proposed position overlaps with any other particle in 
the simulation box then the move will be automatically rejected. Further con¬ 
straints can be applied, for example, on the coordination number and/or nearest 
neighbor distances and angles 31, 32, isl. Hdl. [ 2 ^ . 


The different approaches that we present here, are based on two different 
ways of calculating the total scattering structure factor , Scale, from the particle 
coordinates. 


2.1. Method I: the ’g(r) route’ (RMC++) 


This approach is based on the one-dimensional Fourier transformation of the 
radial distribution function (RDF). For one component systems, the RDF can 
simply be calculated from the atomic positions as 


gir) 


n{r) 

AVp’ 


( 3 ) 


where n(r) is the number of atoms at a distance between r and r + Ar from 
a central atom, AV is the volume of a spherical shell between r and r + Ar and 
p is the number density of the system. 

Liquids and amorphous materials can be considered isotropic beyond nearest- 
neighbor distances so that for switching between the real and reciprocal space, 
a one-dimensional Fourier transform is widely used. Radial distribution func¬ 
tions can be Fourier transformed and weighted for the actual experiment thus 
providing the total scattering structure factor, S{Q). For neutron scattering 
measurements and one component systems, the appropriate Fourier transform 
is given by 


S{Q) = 1 + 


ATrp{b)‘‘ 


r[g(r) — l]sin{Qr)dr ^ 


( 4 ) 


where p denotes the number density of the sample, (5) is the neutron scat¬ 
tering length of the atom type in question, Q are the moduli of the reciprocal 
lattice vectors and the integral runs over atomic distances r. In practice, a dis¬ 
crete integration using the so called rectangular method is performed with 
a summation whose upper limit is restricted by the half-length of the simula¬ 
tion box. This method is implemented in, for instance, the RMC+-I- 32, [H, 
RMC.POT[i3 and RMCProfiledl, [ 2 ^ software packages. 


2.2. Method II: the ’crystallography route’ (RMCPOW) 


In contrast to Method I, the ’crystallography route’, implemented by the 
software RMCPOw[i3, is based on the super-cell approximation, repeating 
the ’unit cell’ (i.e., in our case, the simulation cell) in each direction. The 
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total scattering structure factor, S{Q), is calculated using a three-dimensional 
Fourier transformation to the reciprocal space from atomic coordinates. In 
this way, RMCPOW can deal with ordered and disordered systems because 
diffuse (local disorder) and Bragg scattering (crystalline, long range order) are 
both considered. Diffuse intensities, that are assumed to vary smoothly, are 
locally averaged whereas for Bragg intensities the same summation is performed 
without averaging (see Ref.[l^ for details). 

In the case of neutron diffraction, the orientationally averaged structure 
factor [isjl is 

^ Nv<b >2 ^ (S) 

Where N and V are, respectively, the number of atoms and the volume of 
system, Q' are the allowed vectors in the reciprocal cell, and (b) is the average of 
the coherent scattering lengths. The 1/Q' factor stems from the angular inte¬ 
gration over all the possible Q' orientations[l5|. F(Q) contains the correlations 
between scattering nuclei and is given by 

N 

exp(iQRj), (6) 

1=1 

where Rj denotes the position of atom j in the unit cell. 

It is important to point out that no Fourier transformation is involved in 
this scheme and therefore, the usual numerical problems (truncation, aliasing) 
in conjunction with that do not occur. Another thing to notice is that if Eql5] 
was calculated for an isotropic system without periodical long range ordering 
then the vectors could be substituted by their magnitudes and the summation 
could be replaced by integration; that is, eventually, EqlHwould be reproduced. 

3. Calculations performed 

As mentioned before, approaches I and II are tested here on disordered one 
component systems of increasing complexity, from liquid argon to amorphous 
silicon. The differences in terms of structural order can be clearly seen in FiglU 
where the experimental and simulated structure factors for the four systems are 
shown. Note that as the complexity of the test systems increases, new features 
of the ’diffuse’ scattering appear but not any Bragg peak and therefore method 
I (the ’g(r) route’) can also be used. In all cases, simulated and experimental 
data are from neutron diffraction ’measurements’. 

For 1-Ar, the simplest case, modelled data from canonical Monte Carlo simu¬ 
lations at 85K have been included. In this way, the target structure is accurately 
known and we have access to the real RDF and ADF to compare with. In the 
canonical Monte Carlo simulation argon atoms are modelled using Lennard- 
Jones interactions and the parameters of the LJ potential were taken from the 
literature [^. 
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Figure 1: Comparison of the experimental structure factors provided by experiments (black 
line), RMCH-+ (green line) and RMCPOW (red line). From top to bottom and from left 
to right a) liquid argon (modelled data, from Monte Carlo simulation, see text for details) 
b) liquid gallium (experiment from Ref. [2^ c) liquid selenium (experiment from Ref. [29|l ') d) 
amorphous silicon (experiment from Ref. ). 


All simulations have been performed using the RMC++ 181 and RMCPOW[l5| 
free software packages with cubic simulation boxes of side length 32A for 1-Ga 
and a-Si and 50A for 1-Ar and 1-Se. In table [T] the experimental density and the 
effective size of the particles for each test case are shown. For simplicity, as the 
aim of this work is to compare the two approaches, we have chosen a uniform 
value for the ’experimental’ standard deviation of all Q-values, cr=0.001. As it 
can be seen in Fig[Tl where target and simulated structure factors are hardly 
distinguishable, this value produces good quality fits. 


Table 1: Experimental density and rcutoff used for the model systems. 



p (atoms/A^) 

^cutoff (A) 

1-Ar 

0.02125 

2.7 

1-Ga 

0.05197 

2.2 

1-Se 

0.0298 

2.0 

a-Si 

0.04846 

2.2 
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A better comparison of configurations provided by the different RMC strate¬ 
gies can be made by computing the radial distribution functions (RDF) and 
some simpie three body correlation functions. Such comparison is able to reveal 
subtle variations of the structure that result from the different ways of calcu¬ 
lating the total scattering structure factor. As defined above, the RDF can be 
determined from EqjHl For three body correlations we calculate the bond angle 
distribution (’angular distribution function’, ADF) that can be defined as the 
integral of the three body correlation function g^^^{ri,r2,cos9) over the first 
coordination shell: 


f{0) = 16^2 


rl3dri3rl^dr23g{ri2)g^^^ {ris, ^ 23 , cos 9 ), 


( 7 ) 


where we chose as the position of the first minimum of the radial distribu¬ 
tion function in each test system. This function gives the distribution of angles 
between pairs of nearest neighbors with respect to a central atom. The neutron 
scattering lengths have been taken from Ref. 341. 


4. Results and discussion 
•Liquid argon 

We start by presenting results for argon using the modelled data from the 
canonical Monte Carlo simulation. As it can be seen in FigI2 the agreement is 
almost perfect for RDF and ADF. Therefore it is clearly shown that for ’perfect’ 
experimental data, and for a system with purely two body interactions, both 
RMC approaches reproduce the target structure to the level of two- and three- 
body correlations. 

•Liquid gallium 

In this case the agreement between the RMC-I--I- and the RMCPOW radial 
distribution functions and bond angle distributions (see Figl3]) is almost per¬ 
fect. Behind the good quality of the match of RDF-s and ADF-s one finds the 
considerably higher experimental density of 1-Ga in comparison with 1-Ar (see 
table[T]) and the wider Q range (up to 16A“^) and better quality of the neutron 
scattering measurement. As a consequence of the higher density for gallium, 
the two approaches yield smooth RDF and ADF simply using a simulation box 
of side 32 A. 

•Liquid selenium 

For 1-Se, when comparing the radial and bond angle distribution functions 
in Figure m the overall good agreement is apparent, although the look of these 
functions is not as nice as it was for liquid gallium. A first glance at the RDF 
shows that the fluid is rather structured at short distance, because of the two 
covalent bonds of the atoms; this feature, however, does not seem to impose 
longer range ordering. The short period oscillations of the g{r), again, are 
probably due to some residual systematic errors of the experimental S'((5)[29j. 
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Figure 2: Comparison of the simulated radial distribution functions (top) and bond angle 
distribution functions (bottom) provided by canonical Monte Carlo simulations (black line), 
RMC++ (green line) and RMCPOW (red line) for liquid argon. 


For systems like liquid selenium, in which the short range g(r) displays 
significant features, fine long range details of S{Q) cannot be neglected. Also, 
the size of the simulation box affects to the accuracy of both approaches. By 
increasing the simulation box up to SOAmore particle distances are included in 
the RDF calculation for method I and more reciprocal vectors are included in 
the evaluation of S{Q) for method II. This, in turn, implies both the use of a 
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Figure 3: Comparison of the radial distribution function (top) and the bond angle distribution 
function (bottom) for liquid gallium. 


large system size that allows a finer sampling of r-space (method I) or Q-space 
(method II) and the inclusion of a rather long Q-iange in the fitting procedure. 

The bond angle distributions are very similar, exhibiting maxima at the 
same angle values. The small difference for the maximum at 60 degree angle is 
not unexpected for a liquid with a relatively complex structure (twofold coor¬ 
dination and chain like structure). Only imposing some additional constrains, 
a reasonable good prediction for the three body correlations from a RMC sim- 
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Figure 4: Comparison of the radial distribution function (top) and the bond angle distribution 
function (bottom) for liquid selenium. 


ulation can be obtained for systems with this level of order. 

•Amorphous silicon 

The RMC++ and RMCPOW radial distribution functions (see FiglS]) agree 
very well; that is, the highest level of ordering among our test systems has 
not posed particular difficulties to either approaches. Interestingly, the main 
maximum of the RDF resulting from the 'g{r) route’ is slightly sharper than its 
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counterpart. Since the total scattering structure factors belonging to RMC++ 
and RMCPOW run together, it is not possible to assess which RDF is the ’real’ 
one: one must accept that both (only very slightly different) g{r)-s are possible 
solutions. It would also be rather hard to consider differences in terms the ADF- 
s significant: it might just be noticed that the unphysical maximum at the 60 
degree angles is slightly stronger for the RMC++ solutions, whereas the ’real’ 
maximum around the tetrahedral angle is very similar for the two approaches. 



Figure 5: Comparison of the radial distribution function (top) and the bond angle distribution 
function (bottom) for amorphous silicon. 
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5. Summary and Conclusions 


Two routes to the structure factor calculation implemented in the RMC++ [32 , 
[l^ and RMCPOWflsjl algorithms have been considered. They have been suc¬ 
cessfully tested and compared on four model systems with different level of 
(dis)order. The agreement is almost perfect for simple (1-Ar) and non-covalent 
non-simple (1-Ga) liquids even for more ordered systems (1-Se and a-Si), only 
minor differences appear in terms of the RDF-s and ADF-s. 

In terms of their relative efficiency, we found that the ’g{r) route’ shows a 
considerably faster convergence, but it is exposed to Fourier truncation errors. 
Furthermore, this approach cannot be extended to deal with materials with 
Bragg scattering. In contrast, the ’crystallography’ route can be used, at the 
expense of computational time, for a wide range of systems, from simple liquids 
to perfect crystals. 

For two-phase systems of our future concern, porous crystalline materials 
with partially filled pores, i.e. in which crystalline and liquid/disordered phases 
are simultaneously present in the sample, the ’crystallography route’ is the only 
approach. Out of the presently available software, the RMCPOW algorithm 


seems to be the most general choice; RMCProfile|16l| may also be applicable of 
experimental data over extremely wide Q-range are available. This conclusion is 
also supported by a previous RMC study that showed that the ’crystallography 
route’ provided an appropriate description of simple adsorbed fluids in zeolites, 
although in that case simulated target structure factors were usedf^. 
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